UFR 2-13 Test Case: Difference between revisions

From KBwiki
Jump to navigation Jump to search
m (Dave.Ellacott moved page SilverP:UFR 2-13 Test Case to UFR 2-13 Test Case)
 
(45 intermediate revisions by 3 users not shown)
Line 1: Line 1:
 
=Fluid-structure interaction in turbulent flow past cylinder/plate configuration  I (First swiveling mode)=
= A fluid-structure interaction benchmark in turbulent flow (FSI-PfS-1a) =
{{UFRHeader
{{UFRHeader
|area=2
|area=2
Line 7: Line 6:
__TOC__
__TOC__


The following part is divided into two different sections: in the first one numerical phasedresolved
== Flows Around Bodies ==
results obtained for the two configurations (full and subset case) are compared. Based
=== Underlying Flow Regime 2-13 ===
on this evaluation one case is chosen for a parameter study. Then, in the second subsection
 
the numerical phased-averaged results chosen are juxtaposed to the experimental ones in order
= Test Case Study =
to verify their quality.
== Description of the geometrical model and the test section ==
In both simulations (subset and full case) the flow is initialized by assuming the entire structure
 
to be undeformable. In this case the shell attached to the backside of the cylinder acts
FSI-PfS-1a consists of a flexible thin structure with a distinct
like a splitter plate attenuating the generation of a von Karman vortex street behind the cylinder.
thickness clamped behind a fixed rigid non-rotating cylinder installed
Nevertheless, quasi-periodic vortex shedding is still observed with a Strouhal number
in a water channel (see Fig. 1). The
of St fixed to 0.175. Owing to different loads on both sides the structure starts to deflect as
cylinder has a diameter D=0.022m. It is positioned in the
soon as it is released. After a short initial phase, in which the amplitudes of the deflections
middle of the experimental test section with <math>H_c \operatorname{=} H/2 \operatorname{=} 0.120\,m</math> <math>(H_c/D \approx 5.45)</math>, whereas the test section denotes a
successively increase, a new quasi-periodic mode of oscillation is reached. In accordance with
central part of the entire water channel (see
the experiment in the numerical simulations the shell deforms in the first swiveling mode as
Fig. 2). Its center is located at <math>L_c \operatorname{=}
visible in Fig. 1.
  0.077\,m</math> <math>(L_c/D \operatorname{=} 3.5)</math> downstream of the inflow
section. The test section has a length of <math>L \operatorname{=} 0.338\,m</math>
<math>(L/D \approx 15.36)</math>, a height of <math>H \operatorname{=} 0.240\,m</math>
<math>(H/D \approx 10.91)</math> and a width <math>W \operatorname{=} 0.180\,m</math>
<math>(W/D \approx 8.18)</math>. The blocking ratio of the channel is
about <math>9.2\%</math>. The gravitational acceleration <math>g</math> points in
x-direction (see Fig.~\ref{fig:rubber_plate_geom}), i.e. in the
experimental setup this section of the water channel is turned 90
degrees. The deformable structure used in the experiment behind the
cylinder has a length <math>l \operatorname{=} 0.060\,m</math> <math>(l/D \approx 2.72)</math>
and a width <math>w \operatorname{=} 0.177\,m</math> <math>(w/D \approx 8.05)</math>.
Therefore, in the experiment there is a small gap of about <math>1.5
  \times 10^{-3}\,m</math> between the side of the deformable structure and
both lateral channel walls.
The thickness of the plate is <math>h \operatorname{=} 0.0021\,m</math> <math>(h/D
  \approx 0.09)</math>. This thickness is an averaged value. The material
is natural rubber and thus it is difficult to produce a perfectly
homogeneous 2 mm plate. The measurements show that the thickness of
the plate is between 0.002 and 0.0022 m. All parameters of the
geometrical configuration of the FSI-PfS-1a benchmark are summarized as follows:
 
[[File:table3.png]]
 
[[File:FSI-PfS-1a_Benchmark_Rubberplate_geometry0001_new.png]]
 
Fig. 1: Geometrical configuration of the FSI-PfS-1a Benchmark.
 
== Description of the water channel ==
 
 
In order to validate numerical FSI simulations based on reliable
experimental data, the special research unit on FSI (Bungartz et al. (2006, 2010))
designed and constructed a water channel (Göttingen type, see
Fig. 2) for corresponding experiments with a
special concern regarding controllable and precise boundary and
working conditions Gomes et al. (2006, 2010, 2011, 2013). The
channel (2.8 m x 1.3 m x 0.5 m, fluid volume
of 0.9 m³) has a rectangular flow path and includes several
rectifiers and straighteners to guarantee a uniform inflow into the
test section. To allow optical flow measurement systems like
Particle-Image Velocimetry, the test section is optically accessible
on three sides. It possesses the same geometry as the test section
described in Fig. 1. The structure is
fixed on the backplate of the test section and additionally fixed in
the front glass plate. With a 24 kW axial pump a water inflow of up to
<math>u_{\text{max}}=6</math> m/s is possible. To prevent asymmetries the
gravity force is aligned with the x-axis in main flow direction.
 
[[File:waterchannel_new.png]]
 
Fig. 2: Sketch of the flow channel (dimensions given in mm).
 
== Flow parameters ==
 
Several preliminary tests were performed to find the best working
conditions in terms of maximum structure displacement, good
reproducibility and measurable structure frequencies within the
turbulent flow regime.
 
[[File:inflow1.png]]
 
Fig. 3: Experimental displacements of the structure extremity versus the inlet velocity.
 
Fig. 3 depicts the measured extrema of the structure displacement versus the
inlet velocity and Fig. 4 gives the
frequency and Strouhal number as a function of the inlet
velocity. These data were achieved by measurements with the laser
distance sensor explained in Section [[UFR_2-13_Test Case#Laser distance sensor|Laser distance sensor]]. The
entire diagrams are the result of a measurement campaign in which the
inflow velocity was consecutively increased from 0 to <math>2.2</math> m/s. At an inflow velocity of  <math>u_{\text{inflow}}=1.385 </math> m/s
the displacement are symmetrical, reasonably large and well
reproducible. Based on the inflow velocity chosen and the cylinder
diameter the Reynolds number of the experiment is equal to
<math>\text{Re}=30,470</math>.
 
[[File:inflow2.png]]
 
Fig. 4: Experimental measurements of the frequency and the corresponding Strouhal number of the FSI phenomenon versus the inlet velocity.
 
Regarding the flow around the front
cylinder, at this inflow velocity the flow is in the sub-critical
regime. That means the boundary layers are still laminar, but
transition to turbulence takes place in the free shear layers evolving
from the separated boundary layers behind the apex of the
cylinder. Except the boundary layers at the section walls the inflow
was found to be nearly uniform (see
Fig. 5). The velocity components
<math>\overline{u} </math> and  <math>\overline{v}</math> are measured with two-component
laser-Doppler velocimetry (LDV) along the y-axis in the middle of the
measuring section at  <math>x/D=4.18 </math> and <math>z/D=0 </math>. It can be assumed that
the velocity component <math>\overline{w} </math>  shows a similar velocity profile
as <math>\overline{v} </math>. Furthermore, a low inflow turbulence level of
<math>\text{Tu}_{\text{inflow}}=\sqrt{\frac{1}{3}~\left(\overline{u'^2}+\overline{v'^2}+\overline{w'^2} \right)}/u_{\text{inflow}}=0.02 </math> is measured. All
experiments were performed with water under standard conditions at
<math>T=20^{\circ}\,C</math>. The flow parameters are summarized in the following table:
 
{| class="wikitable"
|+ Flow parameters
|----
! scope="row" | Inflow velocity
| <math>u_{\text{inflow}}=1.385\,m/s </math>
|----
! scope="row" | Flow density
| <math> \rho_f=1000\,kg/m^3</math>
|----
! scope="row" | Flow dynamic viscosity
| <math> \mu_f=1.0 \times 10^{-3}\,Pa\,s </math>
|}
 
[[File:flow_conditions.png]]
 
Fig. 5 Profiles of the mean streamwise and normal velocity as well as the turbulence level at the inflow section of the water channel.
 
== Material Parameters ==
 
Although the material shows a strong non-linear elastic behavior for
large strains, the application of a linear elastic constitutive law
would be favored, to enable the reproduction of this FSI benchmark by
a variety of different computational analysis codes without the need
of complex material laws. This assumption can be justified by the
observation that in the FSI test case, a formulation for large
deformations but small strains is applicable. Hence, the
identification of the material parameters is done on the basis of the
moderate strain expected and the St. Venant-Kirchhoff constitutive law
is chosen as the simplest hyper-elastic material model.
 
The density of the rubber material can be determined to be
<math>\rho_{\text{rubber plate}}</math>=1360 kg/m<math>^3</math> for a thickness
of the plate h = 0.0021 m. This permits the accurate modeling
of inertia effects of the structure and thus dynamic test cases can be
used to calibrate the material constants. For the chosen material
model, there are only two parameters to be defined: the Young's
modulus E and the Poisson's ratio <math>\nu</math>. In order to avoid
complications in the needed element technology due to
incompressibility, the material was realized to have a Poisson's ratio
which reasonably differs from <math>0.5</math>. Material tests of the
manufacturer and complementary experimental/numerical structure test studies (dynamic and decay test scenarios) indicate that the Young's modulus is E=16 MPa and
the Poisson's ratio is <math>\nu</math>=0.48 (a detailed description of the structure tests is available in De Nayer et al., 2014).
 
{| class="wikitable"
|+ Structure parameters
|----
! scope="row" | Density
| <math>\rho_{\text{rubber plate}}=1360\,kg/m^3</math>
|----
! scope="row" | Young's modulus
| <math>E=16\,MPa</math>
|----
! scope="row" | Poisson's ratio
| <math>\nu=0.48\,</math>
|}
 
= Measuring Techniques =
 
Experimental FSI investigations need to contain fluid and structure
measurements for a full description of the coupling process. Under
certain conditions, the same technique for both disciplines can be
used. The measurements performed by
Gomes et al. (2006, 2010, 2013) used the same camera system for
the simultaneous acquisition of the velocity fields and the structural
deflections. This procedure works well for FSI cases involving laminar
flows and 2D structure deflections. In the present case the structure
deforms slightly three-dimensional with increased cycle-to-cycle
variations caused by turbulent variations in the flow. The applied
measuring techniques, especially the structural side, have to deal
with those changed conditions especially the formation of
shades. Furthermore, certain spatial and temporal resolutions as well
as low measurement errors are requested. Due to the different
deformation behavior a single camera setup for the structural
measurements like in Gomes et al. (2006, 2010, 2013) used was not
practicable. Therefore, the velocity fields were captured by a 2D
Particle-Image Velocimetry (PIV) setup and the structural deflections
were measured with a laser triangulation technique. Both devices are
presented in the next sections.
 
=== Particle-image velocimetry ===
 
A classic Particle-Image Velocimetry (Adrian, 1991) setup consists of a single camera obtaining two
components of the fluid velocity on a planar surface illuminated by a
laser light sheet. Particles introduced into the fluid are following the
flow and reflecting the light during the passage of the light sheet.
By taking two reflection fields in a short time interval <math>\Delta</math>t, the most-likely displacements of several particle groups on an
equidistant grid are estimated by a cross-correlation technique or a
particle-tracking algorithm. Based on a precise preliminary
calibration, with the displacements obtained and the time interval
<math>\Delta</math>t chosen, the velocity field can be calculated. To prevent
shadows behind the flexible structure a second light sheet was used to
illuminate the opposite side of the test section.
 
The phase-resolved PIV-measurements (PR-PIV) were carried out with a
4 Mega-pixel camera (TSI Powerview 4MP, charge-coupled device (CCD)
chip) and a pulsed dual-head Neodym:YAG laser (Litron NanoPIV 200)
with an energy of 200 mJ per laser pulse. The high energy of
the laser allowed to use silver-coated hollow glass spheres (SHGS)
with an average diameter of <math>d_{\text{avg,SHGS}}</math>=10~µm and
a density of <math>\rho_{SHGS}</math> = 1400 kg m<math>^{-3}</math> as tracer
particles. To prove the following behavior of these particles a
Stokes number Sk=1.08 and a particle sedimentation velocity
<math>u_{\text{SHGS}}=2.18 \times 10^{-5}~\text{m/s}</math> is calculated.
With this Stokes number and a particle sedimentation velocity which is
much lower than the expected velocities in the experiments, it is ensured that
the tracers closely follow the fluid flow. The camera takes 12 bit pictures with
a frequency of about 7.0 Hz and a resolution of 1695 x 1211 px with respect to the rectangular size of the test
section. For one phase-resolved position (The processing of the phase-resolved fluid velocity fields
involving the structure deflections is described in
Section [[UFR_2-13_Test Case#Generation of Phase-resolved Data|Generation of Phase-resolved Data]].) 60 to 80
measurements are taken. Preliminary studies with more and fewer
measurements showed that this number of measurements represent a good
compromise between accuracy and effort. The grid has a size of 150 x 138 cells and was calibrated with an average
factor of 126 <math>\mu</math> m/px}, covering a planar flow field of
x/D = -2.36 to 7.26 and y/D = -3.47 to ~3.47 in the middle of the test section at
z/D = 0. The time between the frame-straddled laser
pulses was set to <math>\Delta</math> t=200 <math>\mu</math>s. Laser and camera were
controlled by a TSI synchronizer (TSI 610035) with 1 ns
resolution.
 
=== Laser distance sensor ===
 
Non-contact structural measurements are often based on laser distance
techniques. In the present benchmark case the flexible structure shows
an oscillating frequency of about 7.1 Hz. With the
requirement to perform more than 100 measurements per period, a
time-resolved system was needed. Therefore, a laser triangulation was
chosen because of the known geometric dependencies, the high data
rates, the small measurement range and the resulting higher accuracy
in comparison with other techniques such as laser phase-shifting or
laser interferometry. The laser triangulation uses a laser beam which
is focused onto the object. A CCD-chip located near the laser output
detects the reflected light on the object surface. If the distance of
the object from the sensor changes, also the angle changes and thus
the position of its image on the CCD-chip. From this change in
position the distance to the object is calculated by simple
trigonometric functions and an internal length calibration adjusted to
the applied measurement range. To study simultaneously more than one
point on the structure, a multiple-point triangulation sensor was
applied (Micro-Epsilon scanControl 2750, see
Fig. 6). This sensor uses a matrix of CCD
chips to detect the displacements on up to 640 points along a laser
line reflected on the surface of the structure with a data rate of 800 profiles per second. The laser line was positioned in a
horizontal (x/D = 3.2, see
Fig. 6(a)) and in a vertical
alignment(z/D = 0, see
Fig. 6(b)) and has an accuracy of
40 µm}. Due to the different refraction indices of air,
glass and water a custom calibration was performed to take the
modified optical behavior of the emitted laser beams into account.
 
[[File:structure_sensors_scancontrolonly0001_new.png]]
 
Fig. 6: Setup and alignment of multiple-point laser sensor on the flexible structure in a) z-direction and b) x-direction.
 
= Numerical Simulation Methodology =
 
The applied numerical method relies on an efficient partitioned
coupling scheme developed for dynamic fluid-structure interaction
problems in turbulent flows (Breuer et al, 2012). The fluid flow is
predicted by an eddy-resolving scheme, i.e., the large-eddy simulation
technique. FSI problems very often encounter instantaneous
non-equilibrium flows with large-scale flow structures such as
separation, reattachment and vortex shedding. For this kind of flows
the LES technique is obviously the best
choice (Breuer, 2002). Based on a semi-implicit scheme the LES code
is coupled with a code especially suited for the prediction of shells
and membranes. Thus an appropriate tool for the time-resolved
prediction of instantaneous turbulent flows around light, thin-walled
structures results. Since all details of this methodology were
recently published in Breuer et al, 2012, in the following only a
brief description is provided.
 
== Computational fluid dynamics (CFD) ==
 
Within a FSI application the computational domain is no longer fixed
but changes in time due to the fluid forces acting on the structure.
This temporally varying domain is taken into account by the Arbitrary
Lagrangian-Eulerian (ALE) formulation expressing the conservation
equations for time-dependent volumes and surfaces. Here the filtered
Navier-Stokes equations for an incompressible fluid are solved.  Owing
to the deformation of the grid, extra fluxes appear in the governing
equations which are consistently determined considering the space conservation law (SCL) (Demirdzic 1988 and 1990, Lesoinne, 1996). The SCL is expressed
by the swept volumes of the corresponding cell faces and assures that
no space is lost during the movement of the grid lines. For this
purpose the in-house code FASTEST-3D (Durst et al, 1996a, b)
relying on a three-dimensional finite-volume scheme is used. The
discretization is done on a curvilinear, block-structured body-fitted
grid with collocated variable arrangement. A midpoint rule
approximation of second-order accuracy is used for the discretization
of the surface and volume integrals. Furthermore, the flow variables
are linearly interpolated to the cell faces leading to a second-order
accurate central scheme. In order to ensure the coupling of pressure
and velocity fields on non-staggered grids, the momentum interpolation
technique of Rhie and Chow (1983) is used.
 
A predictor-corrector scheme (projection method) of second-order
accuracy forms the kernel of the fluid solver. In the predictor step
an explicit three substep low-storage Runge-Kutta scheme advances
the momentum equation in time leading to intermediate
velocities. These velocities do not satisfy mass conservation. Thus,
in the following corrector step the mass conservation equation has to
be fulfilled by solving a Poisson equation for the
pressure-correction based on the incomplete LU decomposition method
of Stone (1968). The corrector step is repeated (about 3 to 8
iterations) until a predefined convergence criterion (<math>\Delta{m} <
{O}(10^{-9})</math>) is reached and the final velocities and the
pressure of the new time step are obtained. In Breuer et al (2012) it
is explained that the original pressure-correction scheme applied
for fixed grids has not to be changed concerning the mass conservation
equation in the context of moving grids. Solely in the momentum
equation the grid fluxes have to be taken into account as described
above.
 
In LES the large scales in the turbulent flow field are resolved
directly, whereas the non-resolvable small scales have to be taken
into account by a subgrid-scale model. Here the well-known and most
often used eddy-viscosity model, i.e., the Ssmagorinsky (1963) model
is applied. The filter width is directly coupled to the volume of the
computational cell and a Van Driest damping function ensures a
reduction of the subgrid length near solid walls.  Owing to minor
influences of the subgrid-scale model at the moderate Reynolds number
considered in this study, a dynamic procedure to determine the
Smagorinsky parameter as suggested by Germano et al (1991) was omitted and
instead a well established standard constant <math>C_s = 0.1</math> is used.
 
The CFD prediction determines the forces on the structure and delivers
them to the CSD calculation. In the other direction the CSD prediction
determines displacements at the moving boundaries of the computational
domain for the fluid flow. The task is to adapt the grid of the inner
computational domain based on these displacements at the
interface. For moderate deformations algebraic methods are found to be
a good compromise since they are extremely fast and deliver reasonable
grid point distributions maintaining the required high grid
quality. Thus, the grid adjustment is performed based on a transfinite
interpolation (Thompson et al., 1985). It consists of three shear
transformations plus a tensor-product transformation.
 
== Computational structural dynamics (CSD) ==


= Comparison of numerical results =
The dynamic equilibrium of the structure is described by the momentum
equation given in a Lagrangian frame of reference. Large deformations,
where geometrical non-linearities are not negligible, are allowed
(Hojjat et al, 2010). According  to preliminary structure considerations, a total Lagrangian
formulation in terms of the second Piola-Kirchhoff stress tensor and
the Green-Lagrange strain tensor which are linked by the
St. Venant-Kirchhoff material law is used in the present study.
For the solution of the governing equation the finite-element solver
Carat++, which was developed with an emphasis on the prediction of
shell or membrane behavior, is applied. Carat++ is based on several
finite-element types and advanced solution strategies for form finding
and non-linear dynamic
Problems (Wüchner et al, 2005; Bletzinger et al, 2005; Dieringer et al, 2012). For the
dynamic analysis, different time-integration schemes are available,
e.g., the implicit generalized-<math>\alpha</math> method (Chung et al, 1993). In the
modeling of thin-walled structures, membrane or shell elements are
applied for the discretization within the finite-element model. In the
current case, the deformable solid is modeled with a 7-parameter shell
element. Furthermore, special care is given to prevent locking
phenomena by applying the well-known Assumed Natural Strain (ANS) (Hughes and Tezduyar, 1981; Park and Stanley, 1986) and Enhanced Assumed Strain (EAS) methods (Bischoff et al., 2004).


Two numerical setups are used to run the FSI-PfS-1a simulation: the
Both, shell and membrane elements reflect geometrically reduced
full case and the subset case. These configurations differ regarding
structural models with a two-dimensional representation of the
the geometry and the boundary conditions as described in
mid-surface which can describe the three-dimensional physical
Section "Numerical CFD Setup". The subset case represents a
properties by introducing mechanical assumptions for the thickness
simpler model than the full case requiring less CPU-time (one second
direction. Due to this reduced model additional operations are
real-time is predicted in about 170 hours wall-clock with the subset
required to transfer information between the two-dimensional structure
case on 84 processors and in about 310 hours wall-clock with the full
and the three-dimensional fluid model. Thus in the case of shells, the
case on 142 processors) and thus is worth to be considered. The
surface of the interface is found by moving the two-dimensional
question, however, is which influence these modeling assumptions have on
surface of the structure half of the thickness normal to the surface
the numerical results?
on both sides and the closing of the volume (Bletzinger et al., 2006).. On
these two moved surfaces the exchange of data is performed
consistently with respect to the shell theory (Hojjat et al., 2010).


== Full case vs. subset case ==
== Coupling algorithm ==


Both setups are performed with slightly different material
To preserve the advantages of the highly adapted CSD and CFD codes and
characteristics than defined in Section "Material Parameters":
to realize an effective coupling algorithm, a partitioned but
The Young's modulus is set to E=14 MPa, the thickness of
nevertheless strong coupling approach is chosen. Since LES
the plate is equal to h = 0.002 m, the solid density is
typically requires small time steps to resolve the turbulent flow
<math>\rho_\text{rubber plate}</math>=1425 kg m<math>^{-3}</math> and no structural
field, the coupling scheme relies on the explicit
damping is used. The reason is that this comparison was a preliminary
predictor-corrector scheme forming the kernel of the fluid solver.
study carried out prior to the final definition of the test
case. Because of the similitudes of the values used here and those
defined in Section "Material_parameters" and because of the
large CPU-time requested, the comparison of the numerical results is
not repeated with the parameters defined in
Section "Material parameters.


=== Deflection of the structure ===
Based on the velocity and pressure fields from the corrector step, the
fluid forces resulting from the pressure and the viscous shear
stresses at the interface between the fluid and the structure are
computed. These forces are transferred by a grid-to-grid data
interpolation to the CSD code Carat++ using a conservative
interpolation scheme (Farhat et al., 1998) implemented in the coupling
interface CoMA (Gallinger et al, 2009). Using the fluid forces provided via
CoMA, the finite-element code Carat++ determines the stresses in the
structure and the resulting displacements of the structure. This
response of the structure is transferred back to the fluid solver via
CoMA applying a bilinear interpolation which is a consistent scheme
for four-node elements with bilinear shape functions.


At first the predicted deformation of the structure is analyzed. For
For moderate and high density ratios between the fluid and the
this purpose Fig. 1 depicts an
structure, e.g., a flexible structure in water, the added-mass effect
arbitrarily chosen snapshot of the deformed structure for both cases
by the surrounding fluid plays a dominant role. In this situation a
taken from the quasi-periodic oscillation mode. It is observed that
strong coupling scheme taking the tight interaction between the fluid
the shell in the full case deforms more strongly in z-direction than
and the structure into account, is indispensable. In the coupling
in the subset case. This observation can be explained as follows: the
scheme developed in Breuer et al. (2012) this issue is taken into
full setup has a wider structure and the lateral nodes are exposed to
account by a FSI-subiteration loop which works as follows:
less constraints than in the subset case.


[[File:FSI_PfS-1a_num_subset.png]]
A new time step begins with an estimation of the displacement of the
structure. For the estimation a linear extrapolation is applied taking
the displacement values of two former time steps into account.
According to these estimated boundary values, the entire computational
grid has to be adapted as it is done in each FSI-subiteration
loop. Then the predictor-corrector scheme of the next time step is
carried out and the cycle of the FSI-subiteration loop is
entered. After each FSI-subiteration first the FSI convergence is
checked. Convergence is reached if the <math>L_2</math> norm of the displacement
differences between two FSI-subiterations normalized by the <math>L_2</math> norm
of the changes in the displacements between the current and the last
time step drops below a predefined limit, e.g. <math>\varepsilon_{FSI} =
10^{-4}</math> for the present study. Typically, convergence is not reached
within the first step but requires a few FSI-subiterations (5 to
10). Therefore, the procedure has to be continued on the fluid
side. Based on the displacements on the fluid-structure interface,
which are underrelaxated by a constant factor $\omega$ during the
transfer from the CFD to the CSD solver, the inner computational CFD
grid is adjusted. The key point of the coupling procedure suggested in
Breuer et al. (2012) is that subsequently only the corrector step of
the predictor-corrector scheme is carried out again to obtain a new
velocity and pressure field. Thus the clue is that the pressure is
determined in such a manner that the mass conservation is finally
satisfied. Furthermore, this extension of the predictor-corrector
scheme assures that the pressure forces as the most relevant
contribution to the added-mass effect, are successively updated until
dynamic equilibrium is achieved. In conclusion, instabilities due to
the added-mass effect known from loose coupling schemes are avoided
and the explicit character of the time-stepping scheme beneficial for
LES is still maintained.


Fig. 1: Comparison of the structure deformations in y- and z-direction between the full and subset case
The code coupling tool CoMA is based on the
Message-Passing-Interface (MPI) and thus runs in parallel to the
fluid and structure solver. The communication in-between the codes is
performed via standard MPI commands. Since the parallelization in
FASTEST-3D and Carat++ also relies on MPI, a hierarchical
parallelization strategy with different levels of parallelism is
achieved. According to the CPU-time requirements of the different
subtasks, an appropriate number of processors can be assigned to the
fluid and the structure part. Owing to the reduced structural models
on the one side and the fully three-dimensional highly resolved fluid
prediction on the other side, the predominant portion of the CPU-time
is presently required for the CFD part. Additionally, the communication
time between the codes via CoMA and within the CFD solver takes a
non-negligible part of the computational resources.


In order to quantify these displacement variations along the z-axis in
== Numerical CFD Setup ==
the full case, three characteristic points on the structure in three
parallel planes depicted in Fig. 2(c) are
chosen: one plane is set in the middle of the structure, the others
are shifted <math>\pm 60 mm</math> in the spanwise direction. All three points
are not located directly on the shell extremity but at a distance of 9
mm from the extremity. This choice is motivated by the planned
comparison with the measured data (Section Comparison between numerical and experimental results)
and the limitation in the experiment. The laser distance sensor does
not allow to follow the structure extremity and thus points at a
certain distance from the tail are chosen. The dimensionless
y-displacements <math>U_y^* = U_y / D</math> at these three points are
monitored as shown in Fig. 2(a). The following
observation can be made: 1. The displacements are in phase. 2. Local
differences between the curves are observed in the extrema. 3. These
variations are, however, not constant in time. In other words the
displacement in one plane is not always bigger than another. The
variations reflect some kind of waves in the structure that move in
the spanwise direction. Comparing those three raw signals with the
z-averaged displacements depicted in
Fig. 2(b), a maximal difference of 5%
regarding the extrema is noticed. Hence the variations are small. The
corresponding z-variations of the subset case are even smaller (<0.5%). Therefore, it was decided to continue the analysis by averaging
both cases in z-direction.


The next step is to compare the structure deformations obtained with
For the CFD prediction of the flow two different block-structured
the full and the subset case. Fig. 2(b)
grids either for a subset of the entire channel (w/l=1) or for
shows the dimensionless y-displacements of both cases. Notice that by
the full channel but without the gap between the flexible structure
the averaging procedure in z-direction the 3D-problem is reduced to a
and the side walls (w/l=2.95) are used. In the first
2D-problem. The frequencies are identically predicted in both cases
case the entire grid consists of about 13.5 million control volumes
(<math>f_{{FSI}_{\text{num}}} = 6.96</math> Hz and <math>\text{St}_{\text{num}} =
(CVs), whereas 72 equidistant CVs are applied in the spanwise
0.11</math>). Minor differences appear in the extrema of the raw signals
direction. For the full geometry the grid possesses about 22.5 million
presented in Fig. 2(b). As before these
CVs. In this case starting close to both channel walls the grid is
variations are not constant in time and thus the maximal values are
stretched geometrically with a stretching factor 1.1 applying in
found irregularly for either the full or the subset case. As a
total 120 CVs with the first cell center positioned at a distance of
consequence the comparison of the phase-averaged displacement signal
<math>\Delta</math>z/D=1.7 x 10<math>^{-2}</math>.
(see Fig. 2(d)) shows no significant changes
between both cases and the coefficient of determination <math>R^2 = 1 -
\sum_i \left( U_{{y}_i}^* - \hat{U_{{y}_i}^*}\right)^2 / \sum_i \left(U_{{y}_i}^* - \overline{U_y^*}\right)^2</math> of the calculated mean
phase is close to unity (0.9869 for the full case and 0.9782 for the
subset case). <math>\hat{U_{{y}_i}^*}</math> denotes the estimated mean value of
<math>U_y^*</math> for the point i. <math>\overline{U_y^*}</math> is the mean value of all
the displacements. The standard deviation for each point of the
averaged phase is also computed: the maximum for the full case is
0.055 (dimensionless) and for the subset case 0.065
(dimensionless). These values are small compared to the signal, which
is another indication for the reliability of the averaged phase. The
subset case predicts structure deformations very similar to the full
case. In order to check if the FSI results are quasi identical for the
full and the subset case, the phase-resolved flow field has to be
additionally taken into account.


[[File:FSI-PfS-1a_num_full_subset_case.png]]


Fig. 2 Comparison of the structure deformations in y- and z-direction between the full and subset case.
[[File:Benchmark_FSI-PfS-1_full_and_subset_case.jpg]]


=== Phase-resolved flow field ===
Fig. 7: X-Y cross-section of the grid used for the simulation (Note that only every fourth grid line in each direction is displayed here).


The phase-averaging process described in
Section "Generation_of_phase-resolved_data" delivers the
phase-resolved flow fields for the full and the subset case. In order
to compare them just two representative phase-averaged positions of
the FSI problem are chosen to limit this subsection.
Figure 3 shows the flow field in the
vicinity of the shell during its maximal deformation at t=T/4 and Fig. 4 depicts it close to
its undeformed position at t=T, where T denotes the period
time of the phase-averaged signal. The figures display the contours of
the phase-averaged streamwise and transverse velocity
components. Furthermore, the local error of the velocity magnitude
defined by the deviation between the absolute values of the velocity
vector of both cases normalized by the inflow velocity
<math>u_\text{inflow}</math> is depicted. For both positions the results obtained
for the subset and full case are nearly
identical. Figures 3(e)
and 4(e) underline that the local error
of the velocity magnitude between both cases is about zero everywhere
except in the region near the structure. For the position t=T/4 (Fig. 3(e)) small local errors are
located behind the structure in the vortex shedding region. For the
position t=T (Fig. 4(e)) the
phase-averaged position of the shell for the subset case differs
slightly from the one of the full case. Since the flow field is
rapidly changing during the vortex shedding process, this minor
deviation in the phase-angle explains the small local errors observed
near the structure and in the shear layer.


[[File:FSI-PfS-1a_num_flow_1.png]]
The gap between the elastic structure and the walls is not taken into
account in the numerical model and thus the width of the channel is
set to w instead of W. The stretching factors are kept below 1.1 with the first cell
center located at a distance of <math>\Delta</math>y/D=9 x 10<math>^{-4}</math> from
the flexible structure. Based on the wall shear stresses on the
flexible structure the average y<math>^+</math> values are predicted to be below
0.8, mostly even below 0.5. Thus, the viscous sublayer on the
elastic structure and the cylinder is adequately resolved. Since the
boundary layers at the upper and lower channel walls are not
considered, no grid clustering is required here.


Fig. 3 Comparison of the results for the full and subset case; phase-averaged data at t=T/4.
On the CFD side no-slip boundary conditions are applied at the rigid
front cylinder and at the flexible structure. Since the resolution of
the boundary layers at the channel walls would require the bulk of the
CPU-time, the upper and lower channel walls are assumed to be slip
walls. Thus the blocking effect of the walls is maintained without
taking the boundary layers into account. At the inlet a constant
streamwise velocity is set as inflow condition without adding any
perturbations. The choice of zero turbulence level is based on the
consideration that such small perturbations imposed at the inlet will
generally not reach the cylinder due to the coarseness of the grid at
the outer boundaries. Ther


[[File:FSI-PfS-1a_num_flow_2.png]]
== Numerical CSD Setup ==


Fig. 4 Comparison of the results for the full and subset case; phase-averaged data at t=T.
Motivated by the fact that in the case of LES frequently a domain
modeling based on periodic boundary conditions at the lateral walls is
used to reduce the CPU-time requirements, this special approach was
also investigated for the FSI test case. As a consequence, there are two different
structure meshes used: For the CSD prediction of the case with a
subset of the full channel the elastic structure is resolved by the
use of 10 x 10 quadrilateral four-node 7-parameter shell elements. For the
case discretizing the entire channel, 10 quadrilateral four-node shell
elements are used in the main flow direction and 30 in the spanwise
direction.


The comparison of the phase-averaged flow fields shows no significant
On the CSD side, the flexible plate is loaded on the top and bottom
changes between both cases. The subset case predicts the
surface by the fluid forces, which are transferred from the fluid mesh
phase-averaged flow field very similar to the full case. As said
to the structure mesh. These Neumann boundary conditions for the
before, the subset setup is simpler and less expensive in
structure reflect the coupling conditions. Concerning the Dirichlet
CPU-time. Therefore, the subset case is very interesting in order to
boundary conditions, the four edges need appropriate support modeling:
simulate the present test case using LES.
on the upstream side at the rigid cylinder a clamped support is
realized and all degrees of freedom are equal to zero. On the opposite
downstream trailing-edge side, the rubber plate is free to move and
all nodes have the full set of six degrees of freedom. The edges which
are aligned to the main flow direction need different boundary
condition modeling, depending on whether the subset or the full case
is computed:


= Sensitivity study for the subset case =
For the subset case due to the fluid-motivated periodic boundary
conditions, periodicity for the structure is correspondingly assumed
for consistency reasons. As it turns out, this assumption seems to
hold for this specific benchmark configuration and its deformation
pattern which has strong similarity with an oscillation in the first
eigenmode of the plate. Hence, this modeling approach may be used for
the efficient processing of parameter studies, e.g., to evaluate the
sensitivity of the FSI simulations with respect to slight variations
in model parameters shown in a sensitivity study. For
this special type of support modeling, there are always two structure
nodes on the lateral sides (one in a plane z=-w/2 and its twin in
the other plane z=+w/2) which have the same load. These two nodes
must have the same displacements in x- and y-direction and their
rotations have to be identical. Moreover, the periodic boundary
conditions imply that the z-displacement of the nodes on the sides are
forced to be zero.


In order to better understand the test case a comprehensive study on
For the full case the presence of the walls in connection with the
the influence of the three main parameters of the structure (the
small gap implies that there is in fact no constraining effect on the
thickness of the plate h, the density <math>\rho_\text{rubber plate}</math> and
structure, as long as no contact between the plate and the wall takes
the Young's modulus E) was carried out.
place. Out of precise observations in the lab, the possibility of
contact may be disregarded. In principle, this configuration would
lead to free-edge conditions like at the trailing edge. However, the
simulation of the fluid with a moving mesh needs a well-defined mesh
situation at the side walls which made it necessary to tightly connect
the structure mesh to the walls (the detailed representation of the
side edges within the fluid mesh is discarded due to computational
costs and the resulting deformation sensitivity of the mesh in these
regions). Also the displacement in z-direction of the structure nodes
at the lateral boundaries is forced to be zero.


* The thickness of the plate was at first set to h = 0.002 m. However the material is natural rubber and to manufacture a perfectly homogeneous 2 mm plate is not easy. The experimental measurements show that the thickness of the plate varies between 0.002 and 0.0022 m. Therefore, two values of h are tested: the theoretical value of 0.002 m and the average value 0.0021 m.
== Coupling conditions ==


* The density of the plate <math>\rho_\text{rubber plate}</math> is the second parameter. The value of <math>\rho_\text{rubber plate}</math> is determined by a scale and the volume of the structure. Consequently, <math>\rho_\text{rubber plate}</math> also depends on h. With h = 0.002 m <math>\rho_\text{rubber plate}</math> is determined to be equal to 1425 kg m<math>^{-3}</math>. With h = 0.0021 m <math>\rho_\text{rubber plate}</math> is found to be equal to 1360 kg m<math>^{-3}</math>.
For the turbulent flow a time-step size of <math>\Delta t_{f} = 2 \times
10^{-5}s~(\Delta t_{f}^{\ast} = 1.26 \times 10^{-3})</math> in
dimensionless form using <math>u_\text{inflow}</math> and D as reference
quantities) is chosen and the same time-step size is applied for the
structural solver based on the generalized-<math>\alpha</math> method with the
spectral radius <math>\varrho_\infty</math>=1.0, i.e, the Newmark standard
method. For the CFD part this time-step size corresponds to a CFL
number in the order of unity. Furthermore, a constant underrelaxation
factor of <math>\omega</math>=0.5 is considered for the displacements and the
loads are transferred without underrelaxation. In accordance with
previous laminar and turbulent cases in Breuer et al. (2012) the FSI
convergence criterion is set to <math>\varepsilon_{FSI} = 10^{-4}</math> for the
<math>L_2</math> norm of the displacement differences. As estimated from previous
cases Breuer et al. (2012) 5 to 10 FSI-subiterations are required to
reach the convergence criterion.


* The third parameter of the structure is the Young's modulus, because it has an important influence on the modeling of the material. A large spectrum of values for E is tested to evaluate this influence.
After an initial phase in which the coupled system reaches a
statistically steady state, each simulation is carried out for about
4 s real-time corresponding to about 27 swiveling cycles of the
flexible structure.


All the tests were carried out without structural damping and are
For the coupled LES predictions the national supercomputer
summarized in Table 1. The full case
SuperMIG/SuperMUC was used applying either 82 or 140 processors for
used in Section "Full case vs. subset case " and the experimental results are also added as references. Each simulation was
the CFD part of the reduced and full geometry,
done for a time interval of 4 s physical time and comprises about
respectively. Additionally, one processor is required for the coupling
27 swiveling periods. The frequency <math>f_{FSI}</math> of the swiveling mode
code and one processor for the CSD code, respectively.
and the extrema of the mean period of the FSI phenomenon (here the
dimensionless y-displacement <math>U_y^* = U_y / D</math> as explained in
Section "Full case vs. subset case" are compared. Furthermore,
the relative errors between the numerical and experimental values is
given.


[[File:FSI-PfS-1a_parameters_study.png]]
= Generation of Phase-resolved Data =


Tab. 1: Parameter study for the subset case of the FSI test case (without structural damping).
Each flow characteristic of a quasi-periodic FSI problem can be
written as a function <math>f=\bar{f}+\tilde{f}+{f}\prime</math>, where <math>\bar{f}</math>
describes the global mean part, <math>\tilde{f}</math> the quasi-periodic part
and <math>{f}\prime</math> a random turbulence-related
part (Reynolds et al., 1972; Cantwell et al., 1983). This splitting can also be written
in the form <math>f = <f> + f\prime</math>, where <math><f></math> is the
phase-averaged part, i.e., the mean at constant phase. In order to be
able to compare numerical results and experimental measurements, the
irregular turbulent part  <math>f\prime</math> has to be averaged out. This measure is
indispensable owing to the nature of turbulence which solely allows
reasonable comparisons based on statistical data. Therefore, the
present data are phase-averaged to obtain only the phase-resolved
contribution <math><f></math> of the problem, which can be seen as a
representative and thus characteristic signal of the underlying FSI
phenomenon.




The following results and trends can be seen:
The procedure to generate phase-resolved results is the same for the
experiments and the simulations and is also similar to the one
presented in Gomes et al. (2006). The technique can be split up into
three steps:


* By varying the Young's modulus E between 8 and 16 MPa it is possible to control the mode of the FSI phenomenon. Thus E turns out to be the most crucial material parameter. With E smaller than 9 MPa, the system oscillates in the second swiveling mode. With E larger than 12 MPa the structure deflection is dominated by the first bending mode of the structure. For a Young's modulus between 9 and 12 MPa a mode transition phase appears in which both swiveling modes are apparent. In this situation the y-displacements of the plate are no longer quasi-periodic and can not be described by a unique frequency.
* '''Reduce the 3D-problem to a 2D-problem''' - Due to the facts that in the present benchmark the structure deformation in spanwise direction is negligible (see Section [[UFR_2-13_Evaluation#Full case vs. subset case|Full case vs. subset case]] for a discussion on the deviation from 2D) and that the delivered experimental PIV-results are solely available in one x-y-plane, first the 3D-problem is reduced to a 2D-problem. For this purpose the flow field and the plate position in the CFD predictions are averaged in spanwise direction.


* Non-negligible variations in the density (1320 kg m<math>^{-3}</math> <math> \le \rho_\text{rubber plate} \le</math> 1725 kg m<math>^{-3}</math>) for a fixed thickness (h = 0.002 m) and Young's modulus (E = 14 MPa) do not drastically change the results of the frequency and of the mean period extrema. The FSI frequency <math>f_{FSI}</math> slightly decreases with the increase of the density.
* '''Determine n reference positions for the FSI Problem''' - A representative signal of the FSI phenomenon is the history of the y-displacements of the trailing edge of the rubber plate. Therefore, it is used as the trigger signal for this averaging method leading to phase-resolved data. Note that the averaged period of this signal is denoted T. At first, it has to be defined in how many sub-parts the main period of the FSI problem will be divided and so, how many reference positions have to be calculated (for example in the present work n = 23). Then, the margins of each period of the y-displacement curve are determined. In order to do that the intersections between the y-displacement curve and the zero crossings (<math>U_y</math>=0) are looked for and used to limit the periods. Third, each period <math>T_i</math> found is divided into n equidistant sub-parts denoted j.


* Comparing the results for both thicknesses for the range 14 <math> \le E \le </math> 16 MPa, it is obvious that a mild variation of the thickness of the plate (0.1 mm, equivalent to 5%) has a non-negligible influence on the extrema of the mean period and no significant influence on the frequency.
* '''Sort and average the data corresponding to each reference Position''' - The sub-part j of the period <math>T_i</math> corresponds to the sub-part j of the period <math>T_{i+1}</math> and so on. Each data set found in a sub-part j will be averaged with the other sets found in the sub-parts j of all other periods. Finally, data sets of n phase-averaged positions for the representative reference period are achieved.


* Overall the frequency of the FSI phenomenon <math>f_{FSI}</math> is very well predicted (relative error under 2.22%) for all tested parameters leading to the first swiveling mode.
The simulation data containing structure positions, pressure and
velocity fields, are generated every 150 time steps. According to the
frequency observed for the structure and the time-step size chosen
about 50 data sets are obtained per swiveling period. With respect to
the time interval predicted and the number of subparts chosen, the
data for each subpart are averaged from about 50 data sets. A
post-processing program is implemented based on the method described
above. It does not require any special treatment and thus the
aforementioned method to get the phase-resolved results is
straightforward.


* Comparing the results for the density <math>\rho_\text{rubber plate}=1360 kg m^{-3}</math> in the range 14 <math>\le E \le</math> 20 MPa, we observe that the FSI frequency <math>f_{FSI}</math> slightly increases with the Young's modulus and that the displacement extrema decrease.
The current experimental setup consists of the multiple-point triangulation
sensor described in Section [[UFR_2-13_Test Case#Laser distance sensor|Laser distance sensor]] and the
synchronizer of the PIV system. Each measurement pulse of the PIV
system is detected in the data acquisition of the laser distance
sensor, which measures the structure deflection continuously with 800
profiles per second. With this setup, contrary to Gomes et al. (2006),
the periods are not detected during the acquisition but in the
post-processing phase. After the run a specific software based on the
described method mentioned above computes the reference structure
motion period and sorts the PIV data to get the phase-averaged
results. A more detailed description of the phase-averaging method is available in De Nayer et al. (2014).  


In summary, the parameter study shows that the Young's modulus is the
most important parameter: It controls the swiveling mode of the
plate. Furthermore, it can be observed that mild modifications of the
shell thickness have a certain effect on the predicted FSI
phenomenon. Contrarily, this parameter study shows that large
variations of the density do not have major influence on the
predictions. Therefore, errors in the density measurement does not
play an important role. With the support of these extensive
preliminary numerical investigations we can now compare the final
numerical results with the experiment.


----
----
{{ACContribs
{{ACContribs
|authors=G. De Nayer, A. Kalmbach, M. Breuer.
|authors=G. De Nayer, A. Kalmbach, M. Breuer
|organisation= Helmut-Schmidt Universität Hamburg (with support by S. Sicklinger and R. Wüchner from Technische Universität München)
|organisation= Helmut-Schmidt Universität Hamburg (with support by S. Sicklinger and R. Wüchner from Technische Universität München)
}}
}}
Line 230: Line 661:


© copyright ERCOFTAC {{CURRENTYEAR}}
© copyright ERCOFTAC {{CURRENTYEAR}}
<math>Insert formula here</math>

Latest revision as of 12:10, 12 February 2017

Fluid-structure interaction in turbulent flow past cylinder/plate configuration I (First swiveling mode)

Front Page

Description

Test Case Studies

Evaluation

Best Practice Advice

References

Flows Around Bodies

Underlying Flow Regime 2-13

Test Case Study

Description of the geometrical model and the test section

FSI-PfS-1a consists of a flexible thin structure with a distinct thickness clamped behind a fixed rigid non-rotating cylinder installed in a water channel (see Fig. 1). The cylinder has a diameter D=0.022m. It is positioned in the middle of the experimental test section with , whereas the test section denotes a central part of the entire water channel (see Fig. 2). Its center is located at downstream of the inflow section. The test section has a length of , a height of and a width . The blocking ratio of the channel is about . The gravitational acceleration points in x-direction (see Fig.~\ref{fig:rubber_plate_geom}), i.e. in the experimental setup this section of the water channel is turned 90 degrees. The deformable structure used in the experiment behind the cylinder has a length and a width . Therefore, in the experiment there is a small gap of about between the side of the deformable structure and both lateral channel walls. The thickness of the plate is . This thickness is an averaged value. The material is natural rubber and thus it is difficult to produce a perfectly homogeneous 2 mm plate. The measurements show that the thickness of the plate is between 0.002 and 0.0022 m. All parameters of the geometrical configuration of the FSI-PfS-1a benchmark are summarized as follows:

Table3.png

FSI-PfS-1a Benchmark Rubberplate geometry0001 new.png

Fig. 1: Geometrical configuration of the FSI-PfS-1a Benchmark.

Description of the water channel

In order to validate numerical FSI simulations based on reliable experimental data, the special research unit on FSI (Bungartz et al. (2006, 2010)) designed and constructed a water channel (Göttingen type, see Fig. 2) for corresponding experiments with a special concern regarding controllable and precise boundary and working conditions Gomes et al. (2006, 2010, 2011, 2013). The channel (2.8 m x 1.3 m x 0.5 m, fluid volume of 0.9 m³) has a rectangular flow path and includes several rectifiers and straighteners to guarantee a uniform inflow into the test section. To allow optical flow measurement systems like Particle-Image Velocimetry, the test section is optically accessible on three sides. It possesses the same geometry as the test section described in Fig. 1. The structure is fixed on the backplate of the test section and additionally fixed in the front glass plate. With a 24 kW axial pump a water inflow of up to m/s is possible. To prevent asymmetries the gravity force is aligned with the x-axis in main flow direction.

Waterchannel new.png

Fig. 2: Sketch of the flow channel (dimensions given in mm).

Flow parameters

Several preliminary tests were performed to find the best working conditions in terms of maximum structure displacement, good reproducibility and measurable structure frequencies within the turbulent flow regime.

Inflow1.png

Fig. 3: Experimental displacements of the structure extremity versus the inlet velocity.

Fig. 3 depicts the measured extrema of the structure displacement versus the inlet velocity and Fig. 4 gives the frequency and Strouhal number as a function of the inlet velocity. These data were achieved by measurements with the laser distance sensor explained in Section Laser distance sensor. The entire diagrams are the result of a measurement campaign in which the inflow velocity was consecutively increased from 0 to m/s. At an inflow velocity of m/s the displacement are symmetrical, reasonably large and well reproducible. Based on the inflow velocity chosen and the cylinder diameter the Reynolds number of the experiment is equal to .

Inflow2.png

Fig. 4: Experimental measurements of the frequency and the corresponding Strouhal number of the FSI phenomenon versus the inlet velocity.

Regarding the flow around the front cylinder, at this inflow velocity the flow is in the sub-critical regime. That means the boundary layers are still laminar, but transition to turbulence takes place in the free shear layers evolving from the separated boundary layers behind the apex of the cylinder. Except the boundary layers at the section walls the inflow was found to be nearly uniform (see Fig. 5). The velocity components and are measured with two-component laser-Doppler velocimetry (LDV) along the y-axis in the middle of the measuring section at and . It can be assumed that the velocity component shows a similar velocity profile as . Furthermore, a low inflow turbulence level of is measured. All experiments were performed with water under standard conditions at . The flow parameters are summarized in the following table:

Flow parameters
Inflow velocity
Flow density
Flow dynamic viscosity

Flow conditions.png

Fig. 5 Profiles of the mean streamwise and normal velocity as well as the turbulence level at the inflow section of the water channel.

Material Parameters

Although the material shows a strong non-linear elastic behavior for large strains, the application of a linear elastic constitutive law would be favored, to enable the reproduction of this FSI benchmark by a variety of different computational analysis codes without the need of complex material laws. This assumption can be justified by the observation that in the FSI test case, a formulation for large deformations but small strains is applicable. Hence, the identification of the material parameters is done on the basis of the moderate strain expected and the St. Venant-Kirchhoff constitutive law is chosen as the simplest hyper-elastic material model.

The density of the rubber material can be determined to be =1360 kg/m for a thickness of the plate h = 0.0021 m. This permits the accurate modeling of inertia effects of the structure and thus dynamic test cases can be used to calibrate the material constants. For the chosen material model, there are only two parameters to be defined: the Young's modulus E and the Poisson's ratio . In order to avoid complications in the needed element technology due to incompressibility, the material was realized to have a Poisson's ratio which reasonably differs from . Material tests of the manufacturer and complementary experimental/numerical structure test studies (dynamic and decay test scenarios) indicate that the Young's modulus is E=16 MPa and the Poisson's ratio is =0.48 (a detailed description of the structure tests is available in De Nayer et al., 2014).

Structure parameters
Density
Young's modulus
Poisson's ratio

Measuring Techniques

Experimental FSI investigations need to contain fluid and structure measurements for a full description of the coupling process. Under certain conditions, the same technique for both disciplines can be used. The measurements performed by Gomes et al. (2006, 2010, 2013) used the same camera system for the simultaneous acquisition of the velocity fields and the structural deflections. This procedure works well for FSI cases involving laminar flows and 2D structure deflections. In the present case the structure deforms slightly three-dimensional with increased cycle-to-cycle variations caused by turbulent variations in the flow. The applied measuring techniques, especially the structural side, have to deal with those changed conditions especially the formation of shades. Furthermore, certain spatial and temporal resolutions as well as low measurement errors are requested. Due to the different deformation behavior a single camera setup for the structural measurements like in Gomes et al. (2006, 2010, 2013) used was not practicable. Therefore, the velocity fields were captured by a 2D Particle-Image Velocimetry (PIV) setup and the structural deflections were measured with a laser triangulation technique. Both devices are presented in the next sections.

Particle-image velocimetry

A classic Particle-Image Velocimetry (Adrian, 1991) setup consists of a single camera obtaining two components of the fluid velocity on a planar surface illuminated by a laser light sheet. Particles introduced into the fluid are following the flow and reflecting the light during the passage of the light sheet. By taking two reflection fields in a short time interval t, the most-likely displacements of several particle groups on an equidistant grid are estimated by a cross-correlation technique or a particle-tracking algorithm. Based on a precise preliminary calibration, with the displacements obtained and the time interval t chosen, the velocity field can be calculated. To prevent shadows behind the flexible structure a second light sheet was used to illuminate the opposite side of the test section.

The phase-resolved PIV-measurements (PR-PIV) were carried out with a 4 Mega-pixel camera (TSI Powerview 4MP, charge-coupled device (CCD) chip) and a pulsed dual-head Neodym:YAG laser (Litron NanoPIV 200) with an energy of 200 mJ per laser pulse. The high energy of the laser allowed to use silver-coated hollow glass spheres (SHGS) with an average diameter of =10~µm and a density of = 1400 kg m as tracer particles. To prove the following behavior of these particles a Stokes number Sk=1.08 and a particle sedimentation velocity is calculated. With this Stokes number and a particle sedimentation velocity which is much lower than the expected velocities in the experiments, it is ensured that the tracers closely follow the fluid flow. The camera takes 12 bit pictures with a frequency of about 7.0 Hz and a resolution of 1695 x 1211 px with respect to the rectangular size of the test section. For one phase-resolved position (The processing of the phase-resolved fluid velocity fields involving the structure deflections is described in Section Generation of Phase-resolved Data.) 60 to 80 measurements are taken. Preliminary studies with more and fewer measurements showed that this number of measurements represent a good compromise between accuracy and effort. The grid has a size of 150 x 138 cells and was calibrated with an average factor of 126 m/px}, covering a planar flow field of x/D = -2.36 to 7.26 and y/D = -3.47 to ~3.47 in the middle of the test section at z/D = 0. The time between the frame-straddled laser pulses was set to t=200 s. Laser and camera were controlled by a TSI synchronizer (TSI 610035) with 1 ns resolution.

Laser distance sensor

Non-contact structural measurements are often based on laser distance techniques. In the present benchmark case the flexible structure shows an oscillating frequency of about 7.1 Hz. With the requirement to perform more than 100 measurements per period, a time-resolved system was needed. Therefore, a laser triangulation was chosen because of the known geometric dependencies, the high data rates, the small measurement range and the resulting higher accuracy in comparison with other techniques such as laser phase-shifting or laser interferometry. The laser triangulation uses a laser beam which is focused onto the object. A CCD-chip located near the laser output detects the reflected light on the object surface. If the distance of the object from the sensor changes, also the angle changes and thus the position of its image on the CCD-chip. From this change in position the distance to the object is calculated by simple trigonometric functions and an internal length calibration adjusted to the applied measurement range. To study simultaneously more than one point on the structure, a multiple-point triangulation sensor was applied (Micro-Epsilon scanControl 2750, see Fig. 6). This sensor uses a matrix of CCD chips to detect the displacements on up to 640 points along a laser line reflected on the surface of the structure with a data rate of 800 profiles per second. The laser line was positioned in a horizontal (x/D = 3.2, see Fig. 6(a)) and in a vertical alignment(z/D = 0, see Fig. 6(b)) and has an accuracy of 40 µm}. Due to the different refraction indices of air, glass and water a custom calibration was performed to take the modified optical behavior of the emitted laser beams into account.

Structure sensors scancontrolonly0001 new.png

Fig. 6: Setup and alignment of multiple-point laser sensor on the flexible structure in a) z-direction and b) x-direction.

Numerical Simulation Methodology

The applied numerical method relies on an efficient partitioned coupling scheme developed for dynamic fluid-structure interaction problems in turbulent flows (Breuer et al, 2012). The fluid flow is predicted by an eddy-resolving scheme, i.e., the large-eddy simulation technique. FSI problems very often encounter instantaneous non-equilibrium flows with large-scale flow structures such as separation, reattachment and vortex shedding. For this kind of flows the LES technique is obviously the best choice (Breuer, 2002). Based on a semi-implicit scheme the LES code is coupled with a code especially suited for the prediction of shells and membranes. Thus an appropriate tool for the time-resolved prediction of instantaneous turbulent flows around light, thin-walled structures results. Since all details of this methodology were recently published in Breuer et al, 2012, in the following only a brief description is provided.

Computational fluid dynamics (CFD)

Within a FSI application the computational domain is no longer fixed but changes in time due to the fluid forces acting on the structure. This temporally varying domain is taken into account by the Arbitrary Lagrangian-Eulerian (ALE) formulation expressing the conservation equations for time-dependent volumes and surfaces. Here the filtered Navier-Stokes equations for an incompressible fluid are solved. Owing to the deformation of the grid, extra fluxes appear in the governing equations which are consistently determined considering the space conservation law (SCL) (Demirdzic 1988 and 1990, Lesoinne, 1996). The SCL is expressed by the swept volumes of the corresponding cell faces and assures that no space is lost during the movement of the grid lines. For this purpose the in-house code FASTEST-3D (Durst et al, 1996a, b) relying on a three-dimensional finite-volume scheme is used. The discretization is done on a curvilinear, block-structured body-fitted grid with collocated variable arrangement. A midpoint rule approximation of second-order accuracy is used for the discretization of the surface and volume integrals. Furthermore, the flow variables are linearly interpolated to the cell faces leading to a second-order accurate central scheme. In order to ensure the coupling of pressure and velocity fields on non-staggered grids, the momentum interpolation technique of Rhie and Chow (1983) is used.

A predictor-corrector scheme (projection method) of second-order accuracy forms the kernel of the fluid solver. In the predictor step an explicit three substep low-storage Runge-Kutta scheme advances the momentum equation in time leading to intermediate velocities. These velocities do not satisfy mass conservation. Thus, in the following corrector step the mass conservation equation has to be fulfilled by solving a Poisson equation for the pressure-correction based on the incomplete LU decomposition method of Stone (1968). The corrector step is repeated (about 3 to 8 iterations) until a predefined convergence criterion () is reached and the final velocities and the pressure of the new time step are obtained. In Breuer et al (2012) it is explained that the original pressure-correction scheme applied for fixed grids has not to be changed concerning the mass conservation equation in the context of moving grids. Solely in the momentum equation the grid fluxes have to be taken into account as described above.

In LES the large scales in the turbulent flow field are resolved directly, whereas the non-resolvable small scales have to be taken into account by a subgrid-scale model. Here the well-known and most often used eddy-viscosity model, i.e., the Ssmagorinsky (1963) model is applied. The filter width is directly coupled to the volume of the computational cell and a Van Driest damping function ensures a reduction of the subgrid length near solid walls. Owing to minor influences of the subgrid-scale model at the moderate Reynolds number considered in this study, a dynamic procedure to determine the Smagorinsky parameter as suggested by Germano et al (1991) was omitted and instead a well established standard constant is used.

The CFD prediction determines the forces on the structure and delivers them to the CSD calculation. In the other direction the CSD prediction determines displacements at the moving boundaries of the computational domain for the fluid flow. The task is to adapt the grid of the inner computational domain based on these displacements at the interface. For moderate deformations algebraic methods are found to be a good compromise since they are extremely fast and deliver reasonable grid point distributions maintaining the required high grid quality. Thus, the grid adjustment is performed based on a transfinite interpolation (Thompson et al., 1985). It consists of three shear transformations plus a tensor-product transformation.

Computational structural dynamics (CSD)

The dynamic equilibrium of the structure is described by the momentum equation given in a Lagrangian frame of reference. Large deformations, where geometrical non-linearities are not negligible, are allowed (Hojjat et al, 2010). According to preliminary structure considerations, a total Lagrangian formulation in terms of the second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor which are linked by the St. Venant-Kirchhoff material law is used in the present study. For the solution of the governing equation the finite-element solver Carat++, which was developed with an emphasis on the prediction of shell or membrane behavior, is applied. Carat++ is based on several finite-element types and advanced solution strategies for form finding and non-linear dynamic Problems (Wüchner et al, 2005; Bletzinger et al, 2005; Dieringer et al, 2012). For the dynamic analysis, different time-integration schemes are available, e.g., the implicit generalized- method (Chung et al, 1993). In the modeling of thin-walled structures, membrane or shell elements are applied for the discretization within the finite-element model. In the current case, the deformable solid is modeled with a 7-parameter shell element. Furthermore, special care is given to prevent locking phenomena by applying the well-known Assumed Natural Strain (ANS) (Hughes and Tezduyar, 1981; Park and Stanley, 1986) and Enhanced Assumed Strain (EAS) methods (Bischoff et al., 2004).

Both, shell and membrane elements reflect geometrically reduced structural models with a two-dimensional representation of the mid-surface which can describe the three-dimensional physical properties by introducing mechanical assumptions for the thickness direction. Due to this reduced model additional operations are required to transfer information between the two-dimensional structure and the three-dimensional fluid model. Thus in the case of shells, the surface of the interface is found by moving the two-dimensional surface of the structure half of the thickness normal to the surface on both sides and the closing of the volume (Bletzinger et al., 2006).. On these two moved surfaces the exchange of data is performed consistently with respect to the shell theory (Hojjat et al., 2010).

Coupling algorithm

To preserve the advantages of the highly adapted CSD and CFD codes and to realize an effective coupling algorithm, a partitioned but nevertheless strong coupling approach is chosen. Since LES typically requires small time steps to resolve the turbulent flow field, the coupling scheme relies on the explicit predictor-corrector scheme forming the kernel of the fluid solver.

Based on the velocity and pressure fields from the corrector step, the fluid forces resulting from the pressure and the viscous shear stresses at the interface between the fluid and the structure are computed. These forces are transferred by a grid-to-grid data interpolation to the CSD code Carat++ using a conservative interpolation scheme (Farhat et al., 1998) implemented in the coupling interface CoMA (Gallinger et al, 2009). Using the fluid forces provided via CoMA, the finite-element code Carat++ determines the stresses in the structure and the resulting displacements of the structure. This response of the structure is transferred back to the fluid solver via CoMA applying a bilinear interpolation which is a consistent scheme for four-node elements with bilinear shape functions.

For moderate and high density ratios between the fluid and the structure, e.g., a flexible structure in water, the added-mass effect by the surrounding fluid plays a dominant role. In this situation a strong coupling scheme taking the tight interaction between the fluid and the structure into account, is indispensable. In the coupling scheme developed in Breuer et al. (2012) this issue is taken into account by a FSI-subiteration loop which works as follows:

A new time step begins with an estimation of the displacement of the structure. For the estimation a linear extrapolation is applied taking the displacement values of two former time steps into account. According to these estimated boundary values, the entire computational grid has to be adapted as it is done in each FSI-subiteration loop. Then the predictor-corrector scheme of the next time step is carried out and the cycle of the FSI-subiteration loop is entered. After each FSI-subiteration first the FSI convergence is checked. Convergence is reached if the norm of the displacement differences between two FSI-subiterations normalized by the norm of the changes in the displacements between the current and the last time step drops below a predefined limit, e.g. for the present study. Typically, convergence is not reached within the first step but requires a few FSI-subiterations (5 to 10). Therefore, the procedure has to be continued on the fluid side. Based on the displacements on the fluid-structure interface, which are underrelaxated by a constant factor $\omega$ during the transfer from the CFD to the CSD solver, the inner computational CFD grid is adjusted. The key point of the coupling procedure suggested in Breuer et al. (2012) is that subsequently only the corrector step of the predictor-corrector scheme is carried out again to obtain a new velocity and pressure field. Thus the clue is that the pressure is determined in such a manner that the mass conservation is finally satisfied. Furthermore, this extension of the predictor-corrector scheme assures that the pressure forces as the most relevant contribution to the added-mass effect, are successively updated until dynamic equilibrium is achieved. In conclusion, instabilities due to the added-mass effect known from loose coupling schemes are avoided and the explicit character of the time-stepping scheme beneficial for LES is still maintained.

The code coupling tool CoMA is based on the Message-Passing-Interface (MPI) and thus runs in parallel to the fluid and structure solver. The communication in-between the codes is performed via standard MPI commands. Since the parallelization in FASTEST-3D and Carat++ also relies on MPI, a hierarchical parallelization strategy with different levels of parallelism is achieved. According to the CPU-time requirements of the different subtasks, an appropriate number of processors can be assigned to the fluid and the structure part. Owing to the reduced structural models on the one side and the fully three-dimensional highly resolved fluid prediction on the other side, the predominant portion of the CPU-time is presently required for the CFD part. Additionally, the communication time between the codes via CoMA and within the CFD solver takes a non-negligible part of the computational resources.

Numerical CFD Setup

For the CFD prediction of the flow two different block-structured grids either for a subset of the entire channel (w/l=1) or for the full channel but without the gap between the flexible structure and the side walls (w/l=2.95) are used. In the first case the entire grid consists of about 13.5 million control volumes (CVs), whereas 72 equidistant CVs are applied in the spanwise direction. For the full geometry the grid possesses about 22.5 million CVs. In this case starting close to both channel walls the grid is stretched geometrically with a stretching factor 1.1 applying in total 120 CVs with the first cell center positioned at a distance of z/D=1.7 x 10.


Benchmark FSI-PfS-1 full and subset case.jpg

Fig. 7: X-Y cross-section of the grid used for the simulation (Note that only every fourth grid line in each direction is displayed here).


The gap between the elastic structure and the walls is not taken into account in the numerical model and thus the width of the channel is set to w instead of W. The stretching factors are kept below 1.1 with the first cell center located at a distance of y/D=9 x 10 from the flexible structure. Based on the wall shear stresses on the flexible structure the average y values are predicted to be below 0.8, mostly even below 0.5. Thus, the viscous sublayer on the elastic structure and the cylinder is adequately resolved. Since the boundary layers at the upper and lower channel walls are not considered, no grid clustering is required here.

On the CFD side no-slip boundary conditions are applied at the rigid front cylinder and at the flexible structure. Since the resolution of the boundary layers at the channel walls would require the bulk of the CPU-time, the upper and lower channel walls are assumed to be slip walls. Thus the blocking effect of the walls is maintained without taking the boundary layers into account. At the inlet a constant streamwise velocity is set as inflow condition without adding any perturbations. The choice of zero turbulence level is based on the consideration that such small perturbations imposed at the inlet will generally not reach the cylinder due to the coarseness of the grid at the outer boundaries. Ther

Numerical CSD Setup

Motivated by the fact that in the case of LES frequently a domain modeling based on periodic boundary conditions at the lateral walls is used to reduce the CPU-time requirements, this special approach was also investigated for the FSI test case. As a consequence, there are two different structure meshes used: For the CSD prediction of the case with a subset of the full channel the elastic structure is resolved by the use of 10 x 10 quadrilateral four-node 7-parameter shell elements. For the case discretizing the entire channel, 10 quadrilateral four-node shell elements are used in the main flow direction and 30 in the spanwise direction.

On the CSD side, the flexible plate is loaded on the top and bottom surface by the fluid forces, which are transferred from the fluid mesh to the structure mesh. These Neumann boundary conditions for the structure reflect the coupling conditions. Concerning the Dirichlet boundary conditions, the four edges need appropriate support modeling: on the upstream side at the rigid cylinder a clamped support is realized and all degrees of freedom are equal to zero. On the opposite downstream trailing-edge side, the rubber plate is free to move and all nodes have the full set of six degrees of freedom. The edges which are aligned to the main flow direction need different boundary condition modeling, depending on whether the subset or the full case is computed:

For the subset case due to the fluid-motivated periodic boundary conditions, periodicity for the structure is correspondingly assumed for consistency reasons. As it turns out, this assumption seems to hold for this specific benchmark configuration and its deformation pattern which has strong similarity with an oscillation in the first eigenmode of the plate. Hence, this modeling approach may be used for the efficient processing of parameter studies, e.g., to evaluate the sensitivity of the FSI simulations with respect to slight variations in model parameters shown in a sensitivity study. For this special type of support modeling, there are always two structure nodes on the lateral sides (one in a plane z=-w/2 and its twin in the other plane z=+w/2) which have the same load. These two nodes must have the same displacements in x- and y-direction and their rotations have to be identical. Moreover, the periodic boundary conditions imply that the z-displacement of the nodes on the sides are forced to be zero.

For the full case the presence of the walls in connection with the small gap implies that there is in fact no constraining effect on the structure, as long as no contact between the plate and the wall takes place. Out of precise observations in the lab, the possibility of contact may be disregarded. In principle, this configuration would lead to free-edge conditions like at the trailing edge. However, the simulation of the fluid with a moving mesh needs a well-defined mesh situation at the side walls which made it necessary to tightly connect the structure mesh to the walls (the detailed representation of the side edges within the fluid mesh is discarded due to computational costs and the resulting deformation sensitivity of the mesh in these regions). Also the displacement in z-direction of the structure nodes at the lateral boundaries is forced to be zero.

Coupling conditions

For the turbulent flow a time-step size of in dimensionless form using and D as reference quantities) is chosen and the same time-step size is applied for the structural solver based on the generalized- method with the spectral radius =1.0, i.e, the Newmark standard method. For the CFD part this time-step size corresponds to a CFL number in the order of unity. Furthermore, a constant underrelaxation factor of =0.5 is considered for the displacements and the loads are transferred without underrelaxation. In accordance with previous laminar and turbulent cases in Breuer et al. (2012) the FSI convergence criterion is set to for the norm of the displacement differences. As estimated from previous cases Breuer et al. (2012) 5 to 10 FSI-subiterations are required to reach the convergence criterion.

After an initial phase in which the coupled system reaches a statistically steady state, each simulation is carried out for about 4 s real-time corresponding to about 27 swiveling cycles of the flexible structure.

For the coupled LES predictions the national supercomputer SuperMIG/SuperMUC was used applying either 82 or 140 processors for the CFD part of the reduced and full geometry, respectively. Additionally, one processor is required for the coupling code and one processor for the CSD code, respectively.

Generation of Phase-resolved Data

Each flow characteristic of a quasi-periodic FSI problem can be written as a function , where describes the global mean part, the quasi-periodic part and a random turbulence-related part (Reynolds et al., 1972; Cantwell et al., 1983). This splitting can also be written in the form , where is the phase-averaged part, i.e., the mean at constant phase. In order to be able to compare numerical results and experimental measurements, the irregular turbulent part has to be averaged out. This measure is indispensable owing to the nature of turbulence which solely allows reasonable comparisons based on statistical data. Therefore, the present data are phase-averaged to obtain only the phase-resolved contribution of the problem, which can be seen as a representative and thus characteristic signal of the underlying FSI phenomenon.


The procedure to generate phase-resolved results is the same for the experiments and the simulations and is also similar to the one presented in Gomes et al. (2006). The technique can be split up into three steps:

  • Reduce the 3D-problem to a 2D-problem - Due to the facts that in the present benchmark the structure deformation in spanwise direction is negligible (see Section Full case vs. subset case for a discussion on the deviation from 2D) and that the delivered experimental PIV-results are solely available in one x-y-plane, first the 3D-problem is reduced to a 2D-problem. For this purpose the flow field and the plate position in the CFD predictions are averaged in spanwise direction.
  • Determine n reference positions for the FSI Problem - A representative signal of the FSI phenomenon is the history of the y-displacements of the trailing edge of the rubber plate. Therefore, it is used as the trigger signal for this averaging method leading to phase-resolved data. Note that the averaged period of this signal is denoted T. At first, it has to be defined in how many sub-parts the main period of the FSI problem will be divided and so, how many reference positions have to be calculated (for example in the present work n = 23). Then, the margins of each period of the y-displacement curve are determined. In order to do that the intersections between the y-displacement curve and the zero crossings (=0) are looked for and used to limit the periods. Third, each period found is divided into n equidistant sub-parts denoted j.
  • Sort and average the data corresponding to each reference Position - The sub-part j of the period corresponds to the sub-part j of the period and so on. Each data set found in a sub-part j will be averaged with the other sets found in the sub-parts j of all other periods. Finally, data sets of n phase-averaged positions for the representative reference period are achieved.

The simulation data containing structure positions, pressure and velocity fields, are generated every 150 time steps. According to the frequency observed for the structure and the time-step size chosen about 50 data sets are obtained per swiveling period. With respect to the time interval predicted and the number of subparts chosen, the data for each subpart are averaged from about 50 data sets. A post-processing program is implemented based on the method described above. It does not require any special treatment and thus the aforementioned method to get the phase-resolved results is straightforward.

The current experimental setup consists of the multiple-point triangulation sensor described in Section Laser distance sensor and the synchronizer of the PIV system. Each measurement pulse of the PIV system is detected in the data acquisition of the laser distance sensor, which measures the structure deflection continuously with 800 profiles per second. With this setup, contrary to Gomes et al. (2006), the periods are not detected during the acquisition but in the post-processing phase. After the run a specific software based on the described method mentioned above computes the reference structure motion period and sorts the PIV data to get the phase-averaged results. A more detailed description of the phase-averaging method is available in De Nayer et al. (2014).




Contributed by: G. De Nayer, A. Kalmbach, M. Breuer — Helmut-Schmidt Universität Hamburg (with support by S. Sicklinger and R. Wüchner from Technische Universität München)


Front Page

Description

Test Case Studies

Evaluation

Best Practice Advice

References


© copyright ERCOFTAC 2024